
rm(list = setdiff(c("run_start", "run_start_total"), ls()))

df <- readRDS("final_data/data_clean.RDS")

# 1) Function for clusterred standard errors ####

prepTex <- function(mod){
  
  if(class(mod) == "lm"){
    
    out_mod <- lmtest::coeftest(mod, vcov = clubSandwich::vcovCR(mod, type = "CR0", cluster = df$id))
    
    coef_names <- names(out_mod[,1])
    est <- out_mod[,1]
    se <- out_mod[,2]
    pval <- out_mod[,4]
    r2 <- summary(mod)$adj.r.squared
    N <- nobs(mod)
    
    gof <- c(r2, N)
    gof_names <- c("Adj. R Squared", "N")
    gof_decimal <- c(T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
  } else if(class(mod) == "polr"){
    
    out_mod <- lmtest::coeftest(mod, vcov = sandwich::vcovCL(mod, cluster = df_exp$respondent))
    
    coef_names <- names(out_mod[,1])
    est <- out_mod[,1]
    se <- out_mod[,2]
    pval <- out_mod[,4]
    aic <- stats::extractAIC(mod)[2]
    deviance <- mod$deviance
    N <- mod$n
    
    gof <- c(aic, deviance, N)
    gof_names <- c("AIC", "Deviance", "N")
    gof_decimal <- c(T,T,F)
    
    texreg::createTexreg(coef.names = coef_names,
                         coef = est,
                         se = se,
                         pvalues = pval,
                         gof = gof,
                         gof.names = gof_names,
                         gof.decimal = gof_decimal)
  }  else {print("Unknown model")}
  
}


# 2) Descriptive evidence ####

m1 <- lm(Q1~ age + gender + origin_language + nationalism + occupation + religion, df)
m4 <- lm(Q4~ age + gender + origin_language + nationalism + occupation + religion, df)
m3 <- lm(Q3~ age + gender + origin_language + nationalism + occupation + religion, df)
m2 <- lm(Q2~ age + gender + origin_language + nationalism + occupation + religion, df)
m5 <- lm(Q5~ age + gender + origin_language + nationalism + occupation + religion, df)

texreg::screenreg(list(m1,m4, m3, m2,m5),
                  custom.model.names = c("Belonging",
                                         "Voting rights",
                                         "Unemployment benefits",
                                         "More Representated",
                                         "Representation by"))


cm1 <- Q1 ~ age + gender + origin_language + nationalism + occupation + religion
cm2 <- Q2 ~ age + gender + origin_language + nationalism + occupation + religion
cm3 <- Q3 ~ age + gender + origin_language + nationalism + occupation + religion
cm4 <- Q4 ~ age + gender + origin_language + nationalism + occupation + religion
cm5 <- Q5 ~ age + gender + origin_language + nationalism + occupation + religion

##### 2.1) Figure A1 - baseline figure for appendix E ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5),
      cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5),
      cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5)) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         statistic = factor(statistic,
                            levels = c("mm"),
                            labels = c("Marginal mean")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95,
                colour = outcome, shape = outcome)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(feature~., scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") +
  scale_y_continuous(limits = c(3,9), breaks = c(3:9), minor_breaks = NULL) +
  scale_color_discrete("Outcome") + 
  scale_shape_discrete("Outcome") + 
  NULL

ggsave("figures/FigureA1.png", width = 21, height = 17, units = "cm")

##### 2.2) Crohnbach's alpha populist attitudes items ####

data <- select(haven::read_spss("nonsharable_data/popbarow13-12-11.sav") %>% 
                 rename('id' = nomem_encr), 
               pop1, pop2, pop3, pop4, pop5, pop6)

alpha_pop <- psych::alpha(data)

print(alpha_pop)

# 3) Belonging to the people ####

##### 3.1) H1: populist citizens are more exclusionary in their conception of the people than non-populist citizens. ####

h1 <- lm(Q1~ pop + left_right + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

###### 3.1.1) Table 4 ####

texreg::htmlreg(list(prepTex(h1)),
                file = "tables/table4.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement"), 
                custom.model.names = c("Model 1"), single.row = T,
                digits = 2)

summary(lm.beta::lm.beta(h1))

###### 3.1.2) Full model table A2 for appendix G ####

texreg::htmlreg(list(prepTex(h1)),
                file = "tables/tableA2.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates national holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim", 
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background",
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"), 
                custom.note = "%stars. Entries are unstandardise linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 1 uses 'Belonging to the people' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 1"), 
                digits = 2)

##### 3.2) H2: Populist citizens are more likely to exclude politicians from their conception of the people than non-populist citizens. ####

###### 3.2.1) Figure 2 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_mean),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_mean) %>% 
        mutate(pop_mean = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_mean = factor(pop_mean,
                           levels = c("Low", "High", "Difference")),
         y_int = ifelse(pop_mean == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_mean, shape = pop_mean)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .7) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/Figure2.png", width = 14, height = 6, units = "cm")

##### 3.3) H3a: Economic attributes have a stronger effect among left-wing populist citizens ####

###### 3.3.1) Figure 3 ####

df$pop_ideoLWP <- factor(df$pop_ideo,
                         levels = c("Non-Populist", "RWP", "LWP"))

df$pop_ideoLWP2 <- factor(df$pop_ideo,
                          levels = c("RWP", "LWP"))

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = "LWP - RWP") %>% 
        rename("pop_ideoLWP" = pop_ideoLWP2)) %>% 
  filter(feature == "occupation" & pop_ideoLWP %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWP,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP, shape = pop_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/Figure3.png", height = 10, width = 17.5, units = "cm")

##### 3.4) H3b: Cultural attributes have a stronger effect among right-wing populist citizens ####

###### 3.4.1) Figure 4 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = "LWP - RWP") %>% 
        rename("pop_ideoLWP" = pop_ideoLWP2)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop_ideoLWP %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWP,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP, shape = pop_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size= .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/Figure4.png", height = 10, width = 17.5, units = "cm")

# 4) Representation ####

##### 4.1) H1: populist citizens are more exclusionary in their conception of the people than non-populist citizens. ####

h1_2 <- lm(Q2~ pop + left_right + age + gender + origin_language + nationalism + occupation + religion +political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)
h1_5 <- lm(Q5~ pop + left_right + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

###### 4.1.1) Table 4 ####

texreg::htmlreg(list(prepTex(h1_2)),
                file = "tables/table5.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement"), 
                custom.model.names = c("Model 2"), single.row = T,
                digits = 2)

texreg::htmlreg(list(prepTex(h1_5)),
                file = "tables/table6.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement"), 
                custom.model.names = c("Model 3"), single.row = T,
                digits = 2)

###### 4.1.2) Full model table A3 for appendix G ####

texreg::htmlreg(list(prepTex(h1_2)),
                file = "tables/tableA3.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim",
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background", 
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"),
                custom.note = "%stars. Entries are unstandardised linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 2 uses perceptions that this 'individual should be represented more' and Model 3 uses 'Desire to be presented by this individual rather than politicans' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 2"), 
                digits = 2)

###### 4.1.3) Full model table A4 for appendix G ####

texreg::htmlreg(list(prepTex(h1_5)),
                file = "tables/tableA4.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim",
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background", 
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"),
                custom.note = "%stars. Entries are unstandardised linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 2 uses perceptions that this 'individual should be represented more' and Model 3 uses 'Desire to be presented by this individual rather than politicans' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 3"), 
                digits = 2)

##### 4.3) H3a: Economic attributes have a stronger effect among left-wing populist citizens - Appendix G ####

###### 4.3.1) Figures A2 for Appendix F####

df$pop_ideoLWP2 <- factor(df$pop_ideo,
                          levels = c("RWP", "LWP"))

rbind(cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP2),
      cregg::cj(df, cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = BY)) %>% 
  filter(feature == "occupation" & pop_ideoLWP2 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP2 = factor(pop_ideoLWP2,
                              levels = c("RWP", "LWP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP2 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP2, shape = pop_ideoLWP2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom",
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA2.png", height = 10, width = 17.5, units = "cm")

###### 4.3.1) Figures A3 for Appendix F####

rbind(cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP2),
      cregg::cj(df, cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = BY)) %>% 
  filter(feature == "occupation" & pop_ideoLWP2 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP2 = factor(pop_ideoLWP2,
                               levels = c("RWP", "LWP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP2 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP2, shape = pop_ideoLWP2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom",
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA3.png", height = 10, width = 17.5, units = "cm")

##### 4.4) H3b: Cultural attributes have a stronger effect among right-wing populist citizens ####

###### 4.4.1) Figures 5 ####

rbind(cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP2),
      cregg::cj(df, cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = BY)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop_ideoLWP2 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP2 = factor(pop_ideoLWP2,
                              levels = c("RWP", "LWP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP2 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP2, shape = pop_ideoLWP2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom",
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/Figure5.png", height = 10, width = 17.5, units = "cm")

###### 4.4.2) Figures 6 ####

rbind(cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP2),
      cregg::cj(df, cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2) %>% 
        mutate(pop_ideoLWP2 = BY)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop_ideoLWP2 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP2 = factor(pop_ideoLWP2,
                               levels = c("RWP", "LWP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWP2 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWP2, shape = pop_ideoLWP2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/Figure6.png", height = 10, width = 17.5, units = "cm")

# 5) Robustness checks ####

##### 5.1) Appendix I3. Separating the effects of populism and host-ideology ####

df$pop_ideoLWP <- relevel(df$pop_ideo6, "LWnP")

df$pop_ideoRWP <- relevel(df$pop_ideo6, "RWnP")

df$pop_ideoLRWP <- relevel(df$pop_ideo6, "LWP")

###### 5.1.1) Figure A19 separate effects of populism and host-ideology for belonging to the people  ####

h3q1 <- lm(Q1~ occupation * pop_ideo6 + age + gender + origin_language + nationalism + religion, df)
h3q4 <- lm(Q4~ occupation * pop_ideo6 + age + gender + origin_language + nationalism + religion, df)
h3q2 <- lm(Q2~ occupation * pop_ideo6 + age + gender + origin_language + nationalism + religion, df)
h3q5 <- lm(Q5~ occupation * pop_ideo6 + age + gender + origin_language + nationalism + religion, df)

texreg::screenreg(list(prepTex(h3q1), prepTex(h3q4), prepTex(h3q2), prepTex(h3q5)),
                  custom.model.names = c("Belonging",
                                         "Voting rights",
                                         "More Representated",
                                         "Represented by"))

h2cm1 <- Q1 ~ age + gender + origin_language + nationalism + occupation + religion
h2cm2 <- Q2 ~ age + gender + origin_language + nationalism + occupation + religion
h2cm4 <- Q4 ~ age + gender + origin_language + nationalism + occupation + religion
h2cm5 <- Q5 ~ age + gender + origin_language + nationalism + occupation + religion

rbind(cregg::cj(df, h2cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h2cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h2cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h2cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h2cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h2cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h2cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h2cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h2cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h2cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h2cm4, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h2cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h2cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h2cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h2cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h2cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h2cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h2cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h2cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h2cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP)) %>%
  filter(feature == "occupation" & pop_ideo6 %in% c("Non-Populist", "CnP", "CP", "RWnP", "RWP", "LWnP", "LWP", "CP - CnP", "LWP - CnP", "LWP - LWnP", "RWP - CnP", "RWP - RWnP", "RWP - LWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower90),
         upper = ifelse(statistic == "mm", upper83, upper90),
         y_int = ifelse(grepl(" - ", pop_ideo6), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  filter(outcome == "Belonging") %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo6, shape = pop_ideo6)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=2, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_manual("Group:", values = 1:20)  + 
  NULL

ggsave("figures/FigureA19.png", width = 24, height = 16, units = "cm")

# 5.1.2) Figure A20 ####

h4q1 <- lm(Q1~ (origin_language + religion) * pop_ideo6 + age + gender + nationalism + occupation, df)
h4q4 <- lm(Q4~ (origin_language + religion) * pop_ideo6 + age + gender + nationalism + occupation, df)
h4q2 <- lm(Q2~ (origin_language + religion) * pop_ideo6 + age + gender + nationalism + occupation, df)
h4q5 <- lm(Q5~ (origin_language + religion) * pop_ideo6 + age + gender + nationalism + occupation, df)

texreg::screenreg(list(prepTex(h4q1), prepTex(h4q4), prepTex(h4q2), prepTex(h4q5)),
                  custom.model.names = c("Belonging",
                                         "Voting rights",
                                         "More Representated",
                                         "Represented by"))

h3cm1 <- Q1 ~ age + gender + origin_language + nationalism + occupation + religion
h3cm2 <- Q2 ~ age + gender + origin_language + nationalism + occupation + religion
h3cm4 <- Q4 ~ age + gender + origin_language + nationalism + occupation + religion
h3cm5 <- Q5 ~ age + gender + origin_language + nationalism + occupation + religion

rbind(cregg::cj(df, h3cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP)) %>%
  filter(feature %in% c("origin_language", "religion") & pop_ideo6 %in% c("Non-Populist", "CnP", "CP", "RWnP", "RWP", "LWnP", "LWP", "CP - CnP", "LWP - CnP", "LWP - LWnP", "RWP - CnP", "RWP - RWnP", "RWP - LWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower90),
         upper = ifelse(statistic == "mm", upper83, upper90),
         y_int = ifelse(grepl(" - ", pop_ideo6), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  filter(outcome == "Belonging") %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo6, shape = pop_ideo6)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=2, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_manual("Group:", values = 1:20)  + 
  NULL

ggsave("figures/FigureA20.png", width = 24, height = 16, units = "cm")


###### 5.2.1) Figures A21 ####

rbind(cregg::cj(df, h3cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP)) %>%
  filter(feature %in% c("origin_language", "religion") & pop_ideo6 %in% c("Non-Populist", "CnP", "CP", "RWnP", "RWP", "LWnP", "LWP", "CP - CnP", "LWP - CnP", "LWP - LWnP", "RWP - CnP", "RWP - RWnP", "RWP - LWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower90),
         upper = ifelse(statistic == "mm", upper83, upper90),
         y_int = ifelse(grepl(" - ", pop_ideo6), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  filter(outcome == "More\nPolitical Representation") %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo6, shape = pop_ideo6)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=2, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_manual("Group:", values = 1:20)  + 
  NULL

ggsave("figures/FigureA21.png", width = 24, height = 16, units = "cm")

# 5.2.2) Figure A22 ####

rbind(cregg::cj(df, h3cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo6),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo6) %>% 
        mutate(pop_ideo6 = BY),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoRWP),
      cregg::cj(df, h3cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLRWP) %>% 
        mutate(pop_ideo6 = BY) %>% 
        select(-pop_ideoLRWP)) %>%
  filter(feature %in% c("origin_language", "religion") & pop_ideo6 %in% c("Non-Populist", "CnP", "CP", "RWnP", "RWP", "LWnP", "LWP", "CP - CnP", "LWP - CnP", "LWP - LWnP", "RWP - CnP", "RWP - RWnP", "RWP - LWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower90),
         upper = ifelse(statistic == "mm", upper83, upper90),
         y_int = ifelse(grepl(" - ", pop_ideo6), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  filter(outcome == "Represented by") %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo6, shape = pop_ideo6)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=2, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_manual("Group:", values = 1:20)  + 
  NULL

ggsave("figures/FigureA22.png", width = 24, height = 16, units = "cm")

##### 5.2) Voting behaviour rather than populist attitudes ####

###### 5.2.1) Belonging to the people ####

####### 5.2.1.1) H1 ####

h1 <- lm(Q1~ pop_vote + left_right + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

summary(h1)

texreg::htmlreg(list(prepTex(h1)),
                file = "tables/tableA6.html",
                custom.coef.map = list("pop_voteYes" = "Populist voting behaviour",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates national holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim", 
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background",
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"), 
                custom.note = "%stars. Entries are unstandardise linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 1 uses 'Belonging to the people' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 1"), 
                digits = 2)

###### 5.2.1.2) H2 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_vote),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_vote) %>% 
        mutate(pop_vote = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_vote = factor(pop_vote,
                           levels = c("Yes", "No", "Difference")),
         y_int = ifelse(pop_vote == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_vote, shape = pop_vote)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populist voting") + 
  scale_shape_discrete("Populist voting") + 
  NULL

ggsave("figures/FigureA12.png", width = 14, height = 6, units = "cm")

###### 5.2.1.3) H3a ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_vote2),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_vote2) %>% 
  mutate(pop_vote2 = "Difference")) %>%
  filter(feature == "occupation") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_vote2 = factor(pop_vote2,
                            levels = c("LWP", "RWP", "Difference"),
                            labels = c("LWP voting", "RWP voting", "Difference")),
         y_int = ifelse(pop_vote2 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(aes(x = forcats::fct_rev(level), y = estimate,
             ymin = lower95, ymax = upper95, 
             colour = pop_vote2, shape = pop_vote2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1)) + 
  coord_flip() +
  facet_grid(outcome ~ statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "Attribute", y = "Estimate") + 
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA13.png", width = 17.5, height = 10, units = "cm")

###### 5.2.1.4) H3b ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_vote2),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_vote2) %>% 
        mutate(pop_vote2 = "Difference")) %>%
  filter(feature %in% c("origin_language", "nationalism", "religion")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_vote2 = factor(pop_vote2,
                            levels = c("LWP", "RWP", "Difference"),
                            labels = c("LWP voting", "RWP voting", "Difference")),
         y_int = ifelse(pop_vote2 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  ggplot(aes(x = forcats::fct_rev(level), y = estimate,
             ymin = lower95, ymax = upper95, 
             colour = pop_vote2, shape = pop_vote2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1)) + 
  coord_flip() +
  facet_grid(outcome ~ statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "Attribute", y = "Estimate") + 
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA14.png", width = 17.5, height = 10, units = "cm")

###### 5.2.2) Representation attitudes ####

###### 5.2.2.1) More representation (H4) ####

rbind(cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_vote2),
      cregg::cj(df, cm2, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_vote2) %>% 
        mutate(pop_vote2 = "Difference")) %>%
  filter(feature %in% c("origin_language", "nationalism", "religion")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_vote2 = factor(pop_vote2,
                            levels = c("LWP", "RWP", "Difference"),
                            labels = c("LWP voting", "RWP voting", "Difference")),
         y_int = ifelse(pop_vote2 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  ggplot(aes(x = forcats::fct_rev(level), y = estimate,
             ymin = lower95, ymax = upper95, 
             colour = pop_vote2, shape = pop_vote2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1)) + 
  coord_flip() +
  facet_grid(outcome ~ statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "Attribute", y = "Estimate") + 
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA15.png", width = 17.5, height = 10, units = "cm")

###### 5.2.2.2) Representation by (H5) ####

rbind(cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_vote2),
      cregg::cj(df, cm5, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_vote2) %>% 
        mutate(pop_vote2 = "Difference")) %>%
  filter(feature %in% c("origin_language", "nationalism", "religion")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_vote2 = factor(pop_vote2,
                            levels = c("LWP", "RWP", "Difference"),
                            labels = c("LWP voting", "RWP voting", "Difference")),
         y_int = ifelse(pop_vote2 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(aes(x = forcats::fct_rev(level), y = estimate,
             ymin = lower95, ymax = upper95, 
             colour = pop_vote2, shape = pop_vote2)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1)) + 
  coord_flip() +
  facet_grid(outcome ~ statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) + 
  labs(x = "Attribute", y = "Estimate") + 
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA16.png", width = 17.5, height = 10, units = "cm")

#### 5.3) Populist attitudes from wave 13 only ####

h1 <- lm(Q1~ popw13 + left_right + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

# 5.3.1) Table A5 ####

texreg::htmlreg(list(prepTex(h1)),
                file = "tables/tableA5.html",
                custom.coef.map = list("popw13" = "Populist attitudes wave 13",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates national holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim", 
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background",
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"), 
                custom.note = "%stars. Entries are unstandardise linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 1 uses 'Belonging to the people' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 1"), 
                digits = 2)

###### 5.3.2) Figure A6 #### 

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_meanw13),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_meanw13) %>% 
        mutate(pop_meanw13 = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_mean = factor(pop_meanw13,
                           levels = c("Low", "High", "Difference")),
         y_int = ifelse(pop_meanw13 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_meanw13, shape = pop_meanw13)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .7) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA6.png", width = 14, height = 6, units = "cm")

###### 5.3.3) Figure A7 ####

df$pop_ideoLWPw13 <- factor(df$pop_ideow13,
                         levels = c("Non-Populist", "RWP", "LWP"))

df$pop_ideoLWP2w13 <- factor(df$pop_ideow13,
                          levels = c("RWP", "LWP"))

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWPw13),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWPw13) %>% 
        mutate(pop_ideoLWPw13 = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2w13) %>% 
        mutate(pop_ideoLWP2w13 = "LWP - RWP") %>% 
        rename("pop_ideoLWPw13" = pop_ideoLWP2w13)) %>% 
  filter(feature == "occupation" & pop_ideoLWPw13 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWPw13,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWPw13 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWPw13, shape = pop_ideoLWPw13)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA7.png", height = 10, width = 17.5, units = "cm")

###### 5.3.5) Figure A8 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWPw13),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWPw13) %>% 
        mutate(pop_ideoLWPw13 = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP2w13) %>% 
        mutate(pop_ideoLWP2w13 = "LWP - RWP") %>% 
        rename("pop_ideoLWPw13" = pop_ideoLWP2w13)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop_ideoLWPw13 %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWPw13,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop_ideoLWPw13 %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideoLWPw13, shape = pop_ideoLWPw13)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size= .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA8.png", height = 10, width = 17.5, units = "cm")

##### 5.4) Anti-immigration sentiments ####

###### 5.4.1) Table A7 ####

h1 <- lm(Q1~ pop + anti_immigration + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

summary(h1)

texreg::htmlreg(list(prepTex(h1)),
                file = "tables/tableA7.html",
                custom.coef.map = list("pop" = "Populist attitudes",
                                       "anti_immigration" = "Anti-immigration sentiments",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates national holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim", 
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background",
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"), 
                custom.note = "%stars. Entries are unstandardise linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 1 uses 'Belonging to the people' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 1"), 
                digits = 2)


###### 5.4.2) Figure A17 ####

df$pop_ideo_imm_PMP <- factor(df$pop_ideo_imm,
                         levels = c("Non-Populist", "AMP", "PMP"))

df$pop_ideo_imm_PMP2 <- factor(df$pop_ideo_imm,
                          levels = c("AMP", "PMP"))

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo_imm_PMP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo_imm_PMP) %>% 
        mutate(pop_ideo_imm_PMP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo_imm_PMP2) %>% 
        mutate(pop_ideo_imm_PMP2 = "PMP - AMP") %>% 
        rename("pop_ideo_imm_PMP" = pop_ideo_imm_PMP2)) %>% 
  filter(feature == "occupation" & pop_ideo_imm_PMP %in% c("AMP", "PMP", "PMP - AMP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideo_imm_PMP = factor(pop_ideo_imm_PMP,
                              levels = c("Non-Populist", "AMP", "PMP", "PMP - NP", "AMP - NP", "PMP - AMP")),
         y_int = ifelse(pop_ideo_imm_PMP %in% c("PMP - NP", "AMP - NP", "PMP - AMP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo_imm_PMP, shape = pop_ideo_imm_PMP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA17.png", height = 10, width = 17.5, units = "cm")

# 5.4.3) Figure A18 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideo_imm_PMP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo_imm_PMP) %>% 
        mutate(pop_ideo_imm_PMP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideo_imm_PMP2) %>% 
        mutate(pop_ideo_imm_PMP2 = "PMP - AMP") %>% 
        rename("pop_ideo_imm_PMP" = pop_ideo_imm_PMP2)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop_ideo_imm_PMP %in% c("AMP", "PMP", "PMP - AMP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideo_imm_PMP = factor(pop_ideo_imm_PMP,
                                   levels = c("Non-Populist", "AMP", "PMP", "PMP - NP", "AMP - NP", "PMP - AMP")),
         y_int = ifelse(pop_ideo_imm_PMP %in% c("PMP - NP", "AMP - NP", "PMP - AMP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_ideo_imm_PMP, shape = pop_ideo_imm_PMP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA18.png", height = 10, width = 17.5, units = "cm")

##### 5.5) Cut-off point populist attitudes at 4 ####

###### 5.5.1) Computing variable ####

df$pop4_ideo <- factor(ifelse(df$left_right_cat3 == "Left" & df$pop_4 == "High", "LWP",
                                ifelse(df$left_right_cat3 == "Right" & df$pop_4 == "High", "RWP",
                                       ifelse(df$left_right_cat3 == "Centre" & df$pop_4 == "High", "CP", "Non-Populist"))),
                         levels = c("Non-Populist", "LWP", "RWP", "CP"))

###### 5.5.2) H2 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_4),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_4) %>% 
        mutate(pop_4 = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_mean = factor(pop_4,
                           levels = c("Low", "High", "Difference")),
         y_int = ifelse(pop_4 == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_4, shape = pop_4)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .7) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA9.png", width = 14, height = 6, units = "cm")

###### 5.5.3) Figure A10 ####

df$pop4_ideoLWP <- factor(df$pop4_ideo,
                            levels = c("Non-Populist", "RWP", "LWP"))

df$pop4_ideoLWP2 <- factor(df$pop4_ideo,
                             levels = c("RWP", "LWP"))

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop4_ideoLWP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop4_ideoLWP) %>% 
        mutate(pop4_ideoLWP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop4_ideoLWP2) %>% 
        mutate(pop4_ideoLWP2 = "LWP - RWP") %>% 
        rename("pop4_ideoLWP" = pop4_ideoLWP2)) %>% 
  filter(feature == "occupation" & pop4_ideoLWP %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop4_ideoLWP,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop4_ideoLWP %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop4_ideoLWP, shape = pop4_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA10.png", height = 10, width = 17.5, units = "cm")

# 5.5.4) Figure A11 ####

rbind(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~pop4_ideoLWP),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop4_ideoLWP) %>% 
        mutate(pop4_ideoLWP = paste0(substr(BY, 1,3), " - NP")),
      cregg::cj(df, cm1, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop4_ideoLWP2) %>% 
        mutate(pop4_ideoLWP2 = "LWP - RWP") %>% 
        rename("pop4_ideoLWP" = pop4_ideoLWP2)) %>% 
  filter(feature %in% c("origin_language", "nationalism", "religion") & pop4_ideoLWP %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop4_ideoLWP,
                              levels = c("Non-Populist", "RWP", "LWP", "LWP - NP", "RWP - NP", "LWP - RWP")),
         y_int = ifelse(pop4_ideoLWP %in% c("LWP - NP", "RWP - NP", "LWP - RWP"), 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Difference in marginal means")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1"),
                          labels = c("Belonging"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop4_ideoLWP, shape = pop4_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  guides(color=guide_legend(nrow=1, byrow=TRUE)) +
  labs(x = "", y = "") +
  scale_color_discrete("Group:") + 
  scale_shape_discrete("Group:") + 
  NULL

ggsave("figures/FigureA11.png", height = 10, width = 17.5, units = "cm")

##### 5.6) Consistency of effects across iterations ####

df$round <- factor(df$round)

###### 5.6.1) Economic dimension ####

bind_rows(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~round),
          cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~round),
          cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~round)) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         statistic = factor(statistic,
                            levels = c("mm"),
                            labels = c("Marginal mean")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q3", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "Unemployment benefits",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  filter(feature == "Occupation") %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95,
                colour = forcats::fct_rev(round), shape = forcats::fct_rev(round))) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(feature~outcome, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") +
  #scale_y_continuous(limits = c(3,9), breaks = c(3:9), minor_breaks = NULL) +
  scale_color_discrete("Round") + 
  scale_shape_discrete("Round") + 
  NULL

ggsave("figures/FigureA4.png", height = 30, width = 25, units = "cm")

###### 5.6.2) Cultural dimension ####

bind_rows(cregg::cj(df, cm1, id = ~ id, estimate = "mm",  h0 = 5, by=~round),
          cregg::cj(df, cm2, id = ~ id, estimate = "mm",  h0 = 5, by=~round),
          cregg::cj(df, cm5, id = ~ id, estimate = "mm",  h0 = 5, by=~round)) %>%
  filter(feature %in% c("origin_language", "nationalism", "religion")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         statistic = factor(statistic,
                            levels = c("mm"),
                            labels = c("Marginal mean")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q4", "Q3", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Voting rights",
                                     "Unemployment benefits",
                                     "More\nPolitical Representation",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95,
                colour = forcats::fct_rev(round), shape = forcats::fct_rev(round))) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(feature~outcome, scales = "free", space = "free") +
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") +
  #scale_y_continuous(limits = c(3,9), breaks = c(3:9), minor_breaks = NULL) +
  scale_color_discrete("Round") + 
  scale_shape_discrete("Round") + 
  NULL

ggsave("figures/FigureA5.png", height = 30, width = 25, units = "cm")

# 6. Appendices ####
##### 6.1) "Right to vote" and " Access to unemployment benefits" ####
###### 6.1.1) Table A11 ####

h1_3 <- lm(Q3~ pop + left_right + age + gender + origin_language + nationalism + occupation + religion +political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)
h1_4 <- lm(Q4~ pop + left_right + age + gender + origin_language + nationalism + occupation + religion + political_interest + political_efficacy + origin_resp + education_resp + age_resp + gender_resp, df)

texreg::htmlreg(list(prepTex(h1_4), prepTex(h1_3)),
                file = "tables/tableA11.html",
                custom.coef.map = list("pop" = "Populism",
                                       "left_right" = "Left-Right Self-Placement",
                                       "age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates national holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim", 
                                       "political_interest" = "Political interest", 
                                       "political_efficacy" = "Political efficacy", 
                                       "origin_respFirst generation western" = "First generation foreign, Western background", 
                                       "origin_respFirst generation non-western" = "First generation foreign, non-Western background", 
                                       "origin_respSecond generation western" = "Second generation foreign, Western background", 
                                       "origin_respSecond generation non-western" = "Second generation foreign, non-Western background",
                                       "education_respIntermediate secondary education" = "Intermediate secondary education", 
                                       "education_respHigher secondary education" = "Hihger secondary education", 
                                       "education_respIntermediate vocational education" = "Intermediate vocational education", 
                                       "education_respHigher vocational training" = "Higher vocational education", 
                                       "education_respUniversity" = "University", 
                                       "age_resp" = "Age respondent", 
                                       "gender_respFemale" = "Gender respondent: Female"),
                custom.note = "%stars. Entries are unstandardise linear regression coefficients. Individual-clustered standard errors in parentheses. All outcome variables are captured on a ten point scale. Model 2 uses perceptions that this 'individual should be represented more' and Model 3 uses 'Desire to be presented by this individual rather than politicans' as an outcome. Reference categories are the following: For Age, the reference is 'Age 25'. For Origin, it is 'Dutch natives'. For occupation, it is 'politicans'. Model 4 excludes this option, the reference category is 'bankers' for this model. For religion, the reference is 'not religious'.",
                custom.model.names = c("Model 4", "Model 5"), 
                digits = 2)

###### 6.1.2) Figure A44 ####

rbind(cregg::cj(df, cm4, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_mean),
      cregg::cj(df, cm4, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_mean) %>% 
        mutate(pop_mean = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_mean = factor(pop_mean,
                           levels = c("Low", "High", "Difference")),
         y_int = ifelse(pop_mean == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q3", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Unemployment", 
                                     "Right to vote",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_mean, shape = pop_mean)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .7) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA44.png", width = 14, height = 6, units = "cm")

# 6.1.3) Figure A45 ####

rbind(cregg::cj(df, cm3, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_mean),
      cregg::cj(df, cm3, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_mean) %>% 
        mutate(pop_mean = "Difference")) %>% 
  filter(level == "Politician") %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         pop_mean = factor(pop_mean,
                           levels = c("Low", "High", "Difference")),
         y_int = ifelse(pop_mean == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal means", "Diffence in marginal means")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q3", "Q4", "Q2", "Q5"),
                          labels = c("Belonging",
                                     "Unemlpoyment\nbenefits", 
                                     "Right to vote",
                                     "More\nRepresentated",
                                     "Represented by"))) %>% 
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95, 
                colour = pop_mean, shape = pop_mean)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower90, ymax = upper90),
                 position = position_dodge(width = 1), linewidth = 1.1) +
  geom_pointrange(position = position_dodge(width = 1), size = .7) + 
  coord_flip() +
  facet_grid(outcome~statistic, scales = "free", space = "free") +
  theme_bw() +
  #theme(legend.position = c(0.85, 0.25)) +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white"), 
        panel.spacing = unit(1, "lines")) +
  labs(x = "", y = "") + 
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA45.png", width = 14, height = 6, units = "cm")


###### 6.1.4) Figure A46 ####

cm13 <- Q1 ~ age + gender + origin_language + nationalism + occupation + religion
cm23 <- Q2 ~ age + gender + origin_language + nationalism + occupation + religion
cm33 <- Q3 ~ age + gender + origin_language + nationalism + occupation + religion
cm43 <- Q4 ~ age + gender + origin_language + nationalism + occupation + religion
cm53 <- Q5 ~ age + gender + origin_language + nationalism + occupation + religion

df$pop_ideoLWP <- factor(df$pop_ideo,
                         levels = c("RWP", "LWP"))

rbind(cregg::cj(df, cm13, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm13, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm23, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm23, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm33, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm33, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm43, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm43, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm53, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm53, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference")) %>% 
  filter(BY %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWP,
                              levels = c("RWP", "LWP", "Difference")),
         y_int = ifelse(pop_ideoLWP == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q2", "Q3", "Q4", "Q5"),
                          labels = c("Belonging",
                                     "More\nRepresentated",
                                     "Unemployment",
                                     "Right to vote",
                                     "Represented by"))) %>% 
  filter(outcome == "Right to vote") %>%
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95,
                colour = pop_ideoLWP, shape = pop_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(feature~statistic, scales = "free", space = "free") +
  labs(x = "", y = "") + 
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white"),
        legend.title = element_blank(),
        legend.key.width = unit(1.5, "cm")) +
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA46.png", height = 30, width = 25, units = "cm")

# 6.1.5) Figure A47 ####

rbind(cregg::cj(df, cm13, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm13, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm23, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm23, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm33, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm33, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm43, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm43, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference"),
      cregg::cj(df, cm53, id = ~ id, estimate = "mm",  h0 = 5, by=~pop_ideoLWP),
      cregg::cj(df, cm53, id = ~ id, estimate = "mm_differences",  h0 = 5, by=~pop_ideoLWP) %>% 
        mutate(pop_ideoLWP = "Difference")) %>% 
  filter(BY %in% c("RWP", "LWP", "LWP - RWP")) %>% 
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error,
         lower = ifelse(statistic == "mm", lower83, lower95),
         upper = ifelse(statistic == "mm", upper83, upper95),
         pop_ideoLWP = factor(pop_ideoLWP,
                              levels = c("RWP", "LWP", "Difference")),
         y_int = ifelse(pop_ideoLWP == "Difference", 0, NA),
         statistic = factor(statistic,
                            levels = c("mm", "mm_difference"),
                            labels = c("Marginal mean", "Diff. marginal mean")),
         feature = factor(feature, 
                          levels = c("origin_language", 
                                     "nationalism", "religion", "occupation", "age", "gender"),
                          labels = c("Origin", "Nat.",
                                     "Religion", "Occupation", "Age", "Gen.")),
         outcome = factor(outcome, 
                          levels = c("Q1", "Q2", "Q3", "Q4", "Q5"),
                          labels = c("Belonging",
                                     "More\nRepresentated",
                                     "Access to social welfare",
                                     "Voting",
                                     "Represented by"))) %>% 
  filter(outcome == "Access to social welfare") %>%
  ggplot(., aes(x=forcats::fct_rev(level), y = estimate,
                ymin = lower95, ymax = upper95,
                colour = pop_ideoLWP, shape = pop_ideoLWP)) +
  geom_hline(aes(yintercept = y_int), lty = "dotted") +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = 1), linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = 1), size = .5) + 
  coord_flip() +
  facet_grid(feature~statistic, scales = "free", space = "free") +
  labs(x = "", y = "") + 
  theme_bw() +
  theme(legend.position = "bottom", 
        strip.background = element_rect(fill = "white"),
        legend.title = element_blank(),
        legend.key.width = unit(1.5, "cm")) +
  scale_color_discrete("Populism") + 
  scale_shape_discrete("Populism") + 
  NULL

ggsave("figures/FigureA47.png", height = 30, width = 25, units = "cm")

#7) Using continuous measures of populism ####

#7.1 Belonging to the people #####

#7.1.1 H2) #####

#7.1.1.1) Figure A23 ####

h2c <- lm(Q1 ~ pop * occupation + age + gender + origin_language + nationalism + religion, df)

marginaleffects::predictions(h2c, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                      occupation = "Politician",
                                                                      age = "45",
                                                                      gender = "female",
                                                                      religion = "Christian",
                                                                      nationalism = "Does not celebrate",
                                                                      origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 linewidth = 1.2) +
  geom_pointrange() + 
  annotate("label",
           x = 4.25, y = 8,
           label = "Profil:\n 45 years old\n Female\n Christian\nNo celebratation\n Politician\nIE-Some dutch",
           fill = "white",
           color = "black",
           size = 3,
           label.size = 0.5,  # border thickness
           label.r = unit(0.15, "lines"))+  # border radius 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Populism", y = "Estimated 'belonging to the people'")

ggsave("figures/FigureA23.png", width = 17.5, height = 10, units = "cm")


#7.1.2 H3) #####

#7.1.2.1) Figure A24 ####

h3c <- lm(Q1 ~ pop * occupation * left_right_cat3 + age + gender + origin_language + nationalism + religion, df)

marginaleffects::predictions(h3c, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                      left_right_cat3 = c("Left", "Right"),
                                                                      occupation = levels(df$occupation),
                                                                      age = "45",
                                                                      gender = "female",
                                                                      religion = "Christian",
                                                                      nationalism = "Does not celebrate",
                                                                      origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(occupation~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'belonging to the people'") +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA24.png", width = 17.5, height = 10, units = "cm")

# Figure A25 ####

marginaleffects::slopes(h3c, variable = "pop", by = c("left_right_cat3", "occupation")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = occupation, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Occupation", y = "Marginal effect of populism on 'belonging to the people'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA25.png", width = 17.5, height = 10, units = "cm")

#7.1.3 H4) #####

#7.1.3.1) Figure A30 ####

h4co <- lm(Q1 ~ pop * left_right_cat3 * (origin_language) + nationalism + religion + age + gender + occupation, df)

marginaleffects::predictions(h4co, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                      left_right_cat3 = c("Left", "Right"),
                                                                      occupation = "Banker",
                                                                      age = "45",
                                                                      gender = "female",
                                                                      religion = "Christian",
                                                                      nationalism = "Does not celebrate",
                                                                      origin_language = levels(df$origin_language))) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(origin_language~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'belonging to the people'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA30.png", width = 17.5, height = 10, units = "cm")

# Figure A31 ####

marginaleffects::slopes(h4co, variable = "pop", by = c("left_right_cat3", "origin_language")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(origin_language), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                                           linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Origin and integration", y = "Marginal effect of populism on 'belonging to the people'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA31.png", width = 17.5, height = 10, units = "cm")

#7.1.3.2) Figure A26 ####

h4cn <- lm(Q1 ~ pop * left_right_cat3 * (nationalism) + origin_language + religion + age + gender + occupation, df)

marginaleffects::predictions(h4cn, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = "Christian",
                                                                       nationalism = levels(df$nationalism),
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(nationalism~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'belonging to the people'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA26.png", width = 17.5, height = 10, units = "cm")

# Figure A27 ####

marginaleffects::slopes(h4cn, variable = "pop", by = c("left_right_cat3", "nationalism")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(nationalism), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Nationalism", y = "Marginal effect of populism on 'belonging to the people'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA27.png", width = 17.5, height = 10, units = "cm")

#7.1.3.3) Figure A28 ####

h4cr <- lm(Q1 ~ pop * left_right_cat3 * (religion) + origin_language + nationalism + age + gender + occupation, df)

marginaleffects::predictions(h4cr, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = levels(df$religion),
                                                                       nationalism = "Does not celebrate",
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(religion~., nrow = 1) + 
  labs(x = "Populism", y = "Estimated 'belonging to the people'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA28.png", width = 17.5, height = 10, units = "cm")

# Figure A29 ####

marginaleffects::slopes(h4cr, variable = "pop", by = c("left_right_cat3", "religion")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(religion), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Religion", y = "Marginal effect of populism on 'belonging to the people'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA29.png", width = 17.5, height = 10, units = "cm")

# 7.2) More political representation ####

#7.2.1) H5 #####

#7.2.1.1) Figure A36 ####

h5co <- lm(Q2 ~ pop * left_right_cat3 * (origin_language) + nationalism + religion + age + gender + occupation, df)

marginaleffects::predictions(h5co, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = "Christian",
                                                                       nationalism = "Does not celebrate",
                                                                       origin_language = levels(df$origin_language))) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(origin_language~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'more political representation'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA36.png", width = 17.5, height = 10, units = "cm")


# Figure A37 ####

marginaleffects::slopes(h5co, variable = "pop", by = c("left_right_cat3", "origin_language")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(origin_language), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Origin and integration", y = "Marginal effect of populism on 'more political representation'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA37.png", width = 17.5, height = 10, units = "cm")

#7.2.1.2) Figure A32 ####

h5cn <- lm(Q2 ~ pop * left_right_cat3 * (nationalism) + origin_language + religion + age + gender + occupation, df)

marginaleffects::predictions(h5cn, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = "Christian",
                                                                       nationalism = levels(df$nationalism),
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(nationalism~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'more political representation'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA32.png", width = 17.5, height = 10, units = "cm")

# Figure A33 ####

marginaleffects::slopes(h5cn, variable = "pop", by = c("left_right_cat3", "nationalism")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(nationalism), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Nationalism", y = "Marginal effect of populism on 'more political representation'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA33.png", width = 17.5, height = 10, units = "cm")

#7.2.1.3) Figure A34 ####

h5cr <- lm(Q2 ~ pop * left_right_cat3 * (religion) + origin_language + nationalism + age + gender + occupation, df)

marginaleffects::predictions(h5cr, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = levels(df$religion),
                                                                       nationalism = "Does not celebrate",
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(religion~., nrow = 1) + 
  labs(x = "Populism", y = "Estimated 'more political representation'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA34.png", width = 17.5, height = 10, units = "cm")

# Figure A35 ####

marginaleffects::slopes(h5cr, variable = "pop", by = c("left_right_cat3", "religion")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(religion), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Religion", y = "Marginal effect of populism on 'belonging to the people'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA35.png", width = 17.5, height = 10, units = "cm")

# 7.3) Representation by ####

#7.3.1) H6 #####

#7.2.1.1) Figure A42 ####

h6co <- lm(Q5 ~ pop * left_right_cat3 * (origin_language) + nationalism + religion + age + gender + occupation, df)

marginaleffects::predictions(h6co, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = "Christian",
                                                                       nationalism = "Does not celebrate",
                                                                       origin_language = levels(df$origin_language))) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(origin_language~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'Representated by'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA42.png", width = 17.5, height = 10, units = "cm")

# Figure A43 ####

marginaleffects::slopes(h6co, variable = "pop", by = c("left_right_cat3", "origin_language")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(origin_language), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Origin and integration", y = "Marginal effect of populism on 'Representated by'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA43.png", width = 17.5, height = 10, units = "cm")

#7.2.1.2) Figure A38 ####

h6cn <- lm(Q5 ~ pop * left_right_cat3 * (nationalism) + origin_language + religion + age + gender + occupation, df)

marginaleffects::predictions(h6cn, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = "Christian",
                                                                       nationalism = levels(df$nationalism),
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(nationalism~., nrow = 2) + 
  labs(x = "Populism", y = "Estimated 'more political representation'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA38.png", width = 17.5, height = 10, units = "cm")

# Figure A39 ####

marginaleffects::slopes(h6cn, variable = "pop", by = c("left_right_cat3", "nationalism")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(nationalism), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Nationalism", y = "Marginal effect of populism on 'Representated by'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA39.png", width = 17.5, height = 10, units = "cm")

#7.2.1.3) Figure A40 ####

h6cr <- lm(Q5 ~ pop * left_right_cat3 * (religion) + origin_language + nationalism + age + gender + occupation, df)

marginaleffects::predictions(h6cr, newdata = marginaleffects::datagrid(pop = c(1:5),
                                                                       left_right_cat3 = c("Left", "Right"),
                                                                       occupation = "Banker",
                                                                       age = "45",
                                                                       gender = "female",
                                                                       religion = levels(df$religion),
                                                                       nationalism = "Does not celebrate",
                                                                       origin_language = "IE-Speaks some Dutch")) %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_line(position = position_dodge(width = .5)) +
  theme_minimal() +
  theme(legend.position = "bottom") + 
  facet_wrap(religion~., nrow = 1) + 
  labs(x = "Populism", y = "Estimated 'Representated by'") +
  scale_color_discrete("Ideology") + 
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA40.png", width = 17.5, height = 10, units = "cm")

# Figure A41 ####

marginaleffects::slopes(h6cr, variable = "pop", by = c("left_right_cat3", "religion")) %>% 
  filter(left_right_cat3 != "Centre") %>% 
  mutate(lower83 = estimate - 1.386 * std.error,
         upper83 = estimate + 1.386 * std.error) %>% 
  ggplot(., aes(x = forcats::fct_rev(religion), y = estimate, ymin = conf.low, ymax = conf.high, colour = left_right_cat3, shape = left_right_cat3)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_linerange(aes(ymin = lower83, ymax = upper83),
                 position = position_dodge(width = .5),
                 linewidth = 1.2) +
  geom_pointrange(position = position_dodge(width = .5)) + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(x = "Religion", y = "Marginal effect of populism on 'Representated by'") +
  coord_flip() +
  scale_colour_discrete("Ideology") +
  scale_shape_discrete("Ideology")

ggsave("figures/FigureA41.png", width = 17.5, height = 10, units = "cm")

# 8) Balance test ####

dep_vars <- c("age_resp", "gender_resp_dummy", "left_right", "pop", "education_resp_dummy")

bal_age <- lm(age_resp ~ age + gender + origin_language + nationalism + occupation + religion, df)
bal_gender <- glm(gender_resp ~ age + gender + origin_language + nationalism + occupation + religion, df, family = "binomial")
bal_left_right <- lm(left_right ~ age + gender + origin_language + nationalism + occupation + religion, df)
bal_pop <- lm(pop ~ age + gender + origin_language + nationalism + occupation + religion, df)
bal_int <- lm(political_interest ~ age + gender + origin_language + nationalism + occupation + religion, df)
bal_eff <- lm(political_efficacy ~ age + gender + origin_language + nationalism + occupation + religion, df)

bal_edu <- nnet::multinom(education_resp ~ age + gender + origin_language + nationalism + occupation + religion, df)
bal_ori <- nnet::multinom(origin_resp ~ age + gender + origin_language + nationalism + occupation + religion, df)

texreg::screenreg(list(bal_age, bal_gender, bal_left_right, bal_pop, bal_int, bal_eff))
texreg::screenreg(list(bal_edu), beside = T)
texreg::screenreg(list(bal_ori), beside = T)

texreg::htmlreg(list(bal_age, bal_gender, bal_left_right, bal_pop, bal_int, bal_eff),
                file = "tables/tableA8.html",
                custom.coef.map = list("age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim"),
                custom.model.names = c("Age", "Gender", "Left-Right", "Populism", "Pol. Interest", "Pol. Efficacy"),
                single.row = F,
                digits = 2)

texreg::htmlreg(list(bal_edu),
                beside = T,
                file = "tables/tableA9.html",
                custom.coef.map = list("age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim"),
                custom.model.names = c("Interm. secondary", "Higher secondary", "Interm. vocational", "Higher vocational", "University"),
                single.row = F,
                digits = 2)

texreg::htmlreg(list(bal_ori),
                beside = T,
                file = "tables/tableA10.html",
                custom.coef.map = list("age45" = "Age 45",
                                       "age65" = "Age 65",
                                       "genderfemale" = "Female",
                                       "origin_languageIE-Speaks some Dutch" = "Irish, speaks some Dutch",
                                       "origin_languageIE-Speaks Dutch well" = "Irish, speaks Dutch well",
                                       "origin_languageMA-Speaks some Dutch" = "Moroccan, speaks some Dutch",
                                       "origin_languageMA-Speaks Dutch well" = "Moroccan, speaks Dutch well",
                                       "origin_languageSY-Speaks some Dutch" = "Syrian, speaks some Dutch",
                                       "origin_languageSY-Speaks Dutch well" = "Syrian, speaks Dutch well",
                                       "nationalismCelebrates holidays" = "Celebrates holidays",
                                       "occupationBanker" = "Banker",
                                       "occupationFarmer" = "Farmer",
                                       "occupationFruitpicker" = "Fruitpicker",
                                       "occupationUniversity lecturer" = "Uni. lecturer",
                                       "occupationMechanic" = "Mechanic",
                                       "occupationCleaner" = "Cleaner",
                                       "occupationUnemployed" = "Unemployed",
                                       "religionChristian" = "Christian",
                                       "religionMuslim" = "Muslim"),
                custom.model.names = c("1st gen. Western", "1st gen. Non-western", "2nd gen. Western", "2nd gen. Non-western"),
                single.row = F,
                digits = 2)
